Código
library(tidyverse)
library(DT)
library(sf)
library(rgdal)
library(raster)
library(terra)
library(leaflet)
library(leaflet.extras)
library(leafem)
library(viridisLite)
library(dplyr)
library(ggplot2)
library(plotly)
library(devtools)En este documento se presenta un análisis de la riqueza de especies de orquídeas en las areas de conservación de Costa Rica. Se utilizan datos geoespaciales y registros de presencia de orquídeas obtenidos de fuentes confiables. Enlaces a las fuentes de datos: - areas de conservación de Costa Rica en Web Feature Service (Sinac): Archivo GeoJSON de areas de conservación de Costa Rica - Registros de presencia de orquídeas de Costa Rica en GBIF: Archivo CSV de registros de presencia de orquídeas de Costa Rica
library(tidyverse)
library(DT)
library(sf)
library(rgdal)
library(raster)
library(terra)
library(leaflet)
library(leaflet.extras)
library(leafem)
library(viridisLite)
library(dplyr)
library(ggplot2)
library(plotly)
library(devtools)areas <-
st_read(
"areas_conservacion_simp_10m.geojson",
quiet = TRUE # para evitar el despliegue de mensajes
)
orquideas <-
st_read(
"orquideas.csv",
options = c(
"X_POSSIBLE_NAMES=decimalLongitude", # columna de longitud decimal
"Y_POSSIBLE_NAMES=decimalLatitude" # columna de latitud decimal
),
quiet = TRUE
)
areas <-
areas |>
st_transform(4326)
st_crs(orquideas) <- 4326orquideas_union_areas <-
st_join(
x = orquideas,
y = dplyr::select(areas, nombre_ac), # selección de columna cod_canton
join = st_within
)
riqueza_especies_orquideas_area <-
orquideas_union_areas |>
st_drop_geometry() |>
group_by(nombre_ac) |>
summarize(riqueza_especies_orquideas = n_distinct(species, na.rm = TRUE))
areas_union_riqueza <-
left_join(
x = areas,
y = dplyr::select(riqueza_especies_orquideas_area, nombre_ac, riqueza_especies_orquideas),
by = "nombre_ac"
) |>
replace_na(list(riqueza_especies_orquideas = 0))
areas_union_riqueza |>
st_drop_geometry() |>
dplyr::select(nombre_ac, riqueza_especies_orquideas) |>
arrange(desc(riqueza_especies_orquideas)) |>
datatable(
colnames = c("Nombre del área de conservación", "Riqueza de especies de orquídeas"),
options = list(
pageLength = 5,
language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
)
)colores_riqueza_especies <-
colorNumeric(
palette = "Reds",
domain = areas_union_riqueza$riqueza_especies_orquideas,
na.color = "transparent"
)
# Paleta de colores de especies
colores_especies <- colorFactor(
palette = viridis(length(unique(orquideas$species))),
domain = orquideas$species
)
# Mapa leaflet
leaflet() |>
setView(
lng = -84.19452,
lat = 9.572735,
zoom = 7) |>
addTiles(group = "Mapa general (OpenStreetMap)") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales (ESRI World Imagery)"
) |>
addPolygons(
data = areas_union_riqueza,
fillColor = ~ colores_riqueza_especies(areas_union_riqueza$riqueza_especies_orquideas),
fillOpacity = 0.8,
color = "black",
stroke = TRUE,
weight = 1.0,
popup = paste(
paste("<strong>Área de conservación:</strong>", areas_union_riqueza$nombre_ac),
paste("<strong>Riqueza de orquídeas:</strong>", areas_union_riqueza$riqueza_especies_orquideas),
sep = '<br/>'
),
group = "Riqueza de orquídeas"
) |>
addScaleBar(
position = "bottomleft",
options = scaleBarOptions(imperial = FALSE)
) |>
addLegend(
position = "bottomleft",
pal = colores_riqueza_especies,
values = areas_union_riqueza$riqueza_especies_orquideas,
group = "Riqueza de especies",
title = "Riqueza de especies"
) |>
addCircleMarkers(
data = orquideas,
stroke = F,
radius = 4,
fillColor = ~colores_especies(orquideas$species),
fillOpacity = 1.0,
popup = paste(
paste0("<strong>Especie: </strong>", orquideas$species),
paste0("<strong>Localidad: </strong>", orquideas$locality),
paste0("<strong>Fecha: </strong>", orquideas$eventDate),
paste0("<strong>Fuente: </strong>", orquideas$institutionCode),
paste0("<a href='", orquideas$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de presencia"
) |>
addLayersControl(
baseGroups = c(
"Mapa general (OpenStreetMap)",
"Imágenes satelitales (ESRI World Imagery)"
),
overlayGroups = c(
"Riqueza de especies",
"Registros de presencia"
)
) |>
addResetMapButton() |>
addSearchOSM() |>
addMouseCoordinates() |>
addFullscreenControl() |>
hideGroup("Registros de presencia") riqueza_especies_orquideas_area <- riqueza_especies_orquideas_area |>
filter(nombre_ac != "")
grafico_barras_ggplot2 <-
riqueza_especies_orquideas_area |>
ggplot(aes(x = reorder(nombre_ac,-riqueza_especies_orquideas), y = riqueza_especies_orquideas)) +
geom_col(
aes(
text = paste0(
"Riqueza de especies de orquídeas en áreas de conservación: ", round(after_stat(y), 2)
)
)
) +
ggtitle("Riqueza de especies de orquídeas en áreas de conservación") +
xlab("Áreas de conservación") +
ylab("Riqueza de especies de orquídeas") +
labs(caption = "Fuente: Sinac") +
theme_dark() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Gráfico de barras plotly
ggplotly(grafico_barras_ggplot2, tooltip = "text") |>
config(locale = 'es')orquideas_top10 <- orquideas_union_areas |>
count(species, sort = TRUE) |>
top_n(10, n)
grafico_barras_ggplot2 <- orquideas_top10 |>
ggplot(aes(x = reorder(species, -n), y = n)) +
geom_bar(stat = "identity", fill = "red",
aes(text = paste0("Cantidad de registros de presencia de especies: ", n))) +
ggtitle("Registros de presencia para las 10 especies de orquídeas con más registros") +
xlab("Especie") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: Sinac") +
theme_dark() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Gráfico de barras plotly
ggplotly(grafico_barras_ggplot2, tooltip = "text") |>
config(locale = 'es')